clc
clear all
close all
% 定义变量范围和名称
nvars = 11;
lb = [270,91,0,0,0,200,0,0,0,0,28];
ub = [667,330,300,400,100,1560,1591,630,1160,18,28];
varNames = {'Cement','water','fly ash','slag','silica-fume','course aggregate','fine aggregate','crumble rubber','chipped rubber','super-plasticizer','curing'};

function [carbon_emission, energy_comsuption,strength] = objectiveFunction(x)
    % 输入 x = [C, W/B, R, S]
    Cement = x(1);         % 水泥用量 (kg/m³)
    water = x(2);       
    flyash = x(3);       
    slag = x(4); 
    sf=x(5);
    ca=x(6);
    fa=x(7);
    crumblerubber=x(8);
    chippedrubber=x(9);
    sp=x(10);
    Curing=x(11);


    % 计算碳排放（假设水和砂的碳排放可忽略）
    carbon_emission = 0.83 * Cement + 0.0003* water+ 0.009 *flyash+0.019*slag+0.024*sf+0.017*ca+0.011*fa+0.0038*crumblerubber+0.0035*chippedrubber+2.388*sp; % 简化的碳排放模型
     energy_comsuption = 4.727* Cement + 0.006* water+ 0.833 *flyash+1.588*slag+1.8*sf+0.022*ca+0.113*fa+0.0135*crumblerubber+0.0135*chippedrubber+18.3*sp
% 强度预测模型（示例：基于经验公式）
    % 假设强度与水泥用量、水胶比、橡胶含量相关
    load('trainedModel.mat')
    predictorNames = trainedModel.RegressionGP.PredictorNames;
    x_table = array2table(x, 'VariableNames', predictorNames);  % 创建输入表格
    x_matrix = table2array(x_table);
    % 使用回归模型预测强度
    strength = trainedModel.predictFcn(x_matrix);  % 强度转为负数，以便NSGA-II最大化; % 示例公式，需替换为实际模型
  
end
function f = multiObjective(x)
    [carbon_emission, energy_comsuption,strength] = objectiveFunction(x);
    f(1) = carbon_emission; % 目标1：最小化碳排放
    f(2) = energy_comsuption;% 目标2：最小化能量消耗
    f(3)=-strength;
end
function [c, ceq] = constraints(x)
    [~,~, strength] = objectiveFunction(x);
    c(1) = 30 - strength; % 强度需 ≥30MPa → 约束 c ≤ 0
   c(2) = strength-40
    ceq=x(1)/3150+x(2)/1000+x(3)/2300+x(4)/2800+x(5)/2300+x(6)/2600+x(7)/2700+x(8)/830+x(9)/1020+x(10)/1300-1;
end
% 设置优化选项
options = optimoptions('gamultiobj', ...
    'PopulationSize', 1000, ...
    'MaxGenerations', 50, ...
    'PlotFcn', {@gaplotpareto}, ...
    'Display', 'iter');

% 运行NSGA-II
[x_optimal, fval] = gamultiobj(@multiObjective, nvars, [], [], [], [], lb, ub, @constraints, options);

% 分析结果
carbon_emissions = fval(:,1);
energy_comsuptions = fval(:,2);
strength=-fval(:,3);
valid_indices = strength >= 30;
valid_solutions = x_optimal(valid_indices, :);
valid_carbon = carbon_emissions(valid_indices);
valid_energy = energy_comsuptions(valid_indices);

% 输出最优解
valid=valid_carbon.^2+valid_energy.^2
[~, min_idx] = min(valid);
optimal_mix = valid_solutions(min_idx, :);
fprintf('optimum mixture:\n');
fprintf('Cement: %.1f kg/m³\n water: %.1f kg/m³\n fly ash: %.1f kg/m³\n slag: %.1f kg/m³\n silica-fume: %.1f kg/m³\n course aggregate: %.1f kg/m³\n  fine aggregate: %.1f kg/m³\n crumble rubber: %.1f kg/m³\n chipped rubber: %.1f kg/m³\n super-plasticizer: %.1f kg/m³\n curing: %.1f days\n', optimal_mix);
fprintf('carbon emmission: %.1f kgCO₂/m³\n compressive strength: %.1f MPa\n', valid_carbon(min_idx), valid_energy(min_idx));

% 绘图
figure;
scatter(carbon_emissions, energy_comsuptions, 'filled');
xlabel('carbon emmission (kgCO₂/m³)');
ylabel('energy comsuption (MJ/m³)');
title('Pareto Front');
hold on;
plot(valid_carbon(min_idx), valid_energy(min_idx), 'ro', 'MarkerSize', 10);

 figure;
scatter3(carbon_emissions, energy_comsuptions, strength, 100, 'filled', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k');
xlabel('Carbon Emissions(MJ/m^3)');
ylabel('Energy Comsuptions(kg CO_2/m^3)');
zlabel('Strength(MPa)');
title('3D Pareto Front');
grid on;